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Abstract 

The gravitational self-force (GSF) and post-Newtonian (PN) schemes are complementary ap¬ 
proximation methods for modelling the dynamics of compact binary systems. Comparison of their 
results in an overlapping domain of validity provides a crucial test for both methods, and can be 
used to enhance their accuracy, e.g. via the determination of previously unknown PN parameters. 
Here, for the first time, we extend such comparisons to noncircular orbits—specifically, to a system 
of two nonspinning objects in a bound (eccentric) orbit. To enable the comparison we use a certain 
orbital-averaged quantity {U) that generalizes Detweiler’s redshift invariant. The functional rela¬ 
tionship {U){Llr,Ll<p), where Hr and are the frequencies of the radial and azimuthal motions, 
is an invariant characteristic of the conservative dynamics. We compute {U){Hr,H^) numerically 
through linear order in the mass ratio q, using a GSF code which is based on a frequency-domain 
treatment of the linearized Einstein equations in the Lorenz gauge. We also derive {U){Hr,H^) 
analytically through 3PN order, for an arbitrary q, using the known near-zone 3PN metric and the 
generalized quasi-Keplerian representation of the motion. We demonstrate that the 0(q) piece of 
the analytical PN prediction is perfectly consistent with the numerical GSF results, and we use 
the latter to estimate yet unknown pieces of the 4PN expression at 0{q). 
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I. INTRODUCTION 


With the Advanced LIGO observatories scheduled to start science runs in 2015 [1], the 
next few years are likely to see hrst direct detections of gravitational waves from astrophysical 
sources. Prime targets are inspiralling and coalescing binary systems of neutron stars and/or 
black holes, with predicted rates that may be as high as a few dozen per observation year [2], 
Theoretical templates of the gravitational waveforms must be developed to enable detection 
and interpretation of the weak signals [3]. The parameter space of these waveforms is too 
large for numerical relativity simulations to cover sufficiently well. Instead, the community 
has been seeking semi-analytical models that can be informed by a judiciously chosen set of 
numerical relativity templates. A leading framework is the effective one-body (EOB) model, 
where the two-body relativistic dynamics is mapped onto a model of (non)geodesic motion 
in an effective spacetime [4-7]. EOB waveforms will play a crucial role in searches based on 
matched filtering, and there is an important need to refine the model, particularly in the 
strong-field regime [8-10]. 

One avenue of refinement is provided by the gravitational self-force (GSF) method, a 
perturbative scheme based on an expansion in the mass ratio of the binary [11-13]. The GSF 
approach is complementary to the post-Newtonian (PN) approximation, a weak-field/small- 
velocity expansion valid for arbitrary mass ratios [14]. Recently, there has been much activity 
in attempt to “synergize” the two schemes. The goal of such cross-cultural studies is three¬ 
fold: to test the two independent approximation schemes—GSF and PN—and help delineate 
their respective domains of validity; to determine yet-unknown high-order expansion terms 
in both approaches (hence improving both approximations); and to help calibrate the EOB 
model across the entire inspiral parameter space. 

To facilitate such studies requires the identification of concrete gauge-invariant physical 
quantities that can be computed using both approaches. A first such quantity was identified 
by Detweiler in 2008 within the GSF framework [15]: the so-called “redshift” variable, de¬ 
fined for strictly circular orbits when dissipation is ignored. The functional relation between 
the redshift and the orbital frequency is a gauge-invariant diagnostic of the conservative 
sector of the binary dynamics. Detweiler made the first successful comparison with the PN 
prediction at 2PN order [15]. This comparison was later extended by Blanchet et al. to 3PN 
order and to even higher orders [16-20]. The calculation of the redshift through linear order 
in the mass ratio was subsequently confirmed by several other GSF computations in different 
gauges [21, 22], which provided an internal consistency check for the GSF formalism. 

Soon after, Barack and Sago considered two more such “conservative” invariant quan¬ 
tities, namely the frequency of the innermost stable circular orbit (ISGO), and the rate of 
periastron advance [23, 24]. These results led to a plethora of comparisons between PN, 
GSF and numerical relativity [25-30] and the subsequent refinement of EOB theory [31-34]. 
More recently, the geodetic spin precession along circular orbits was computed by Dolan et 
al. through linear order in the mass ratio, and the numerical results successfully compared to 
a 3PN-accurate prediction [35]. The results allowed a numerical prediction of the (hitherto 
unknown) 4PN expression for the spin precession. This was later confirmed analytically by 
Bini and Damour [9], who also proceeded to obtain all PN terms up to the 8.5PN order, at 
linear order in the mass ratio. Dolan et al. [36] then presented a computation of the leading 
post-geodesic corrections to certain tidal invariants defined along the orbit. The PN series 
for these tidal invariants were also computed analytically up to 7.5PN order in Ref. [10], 
still at linear order in the mass ratio. 


2 


All synergistic work so far has focused on circular orbits, for simplicity. Here, for the hrst 
time, we extend this program to orbits of arbitrary eccentricity. There are several reasons to 
do so. First, eccentricity provides more “handle” on the strong-held dynamics, giving access 
to new degrees of freedom in the FOB formulation. Second, although most Advanced LIGO 
binaries would have completely circularized by the time they enter the observable frequency 
band, there are scenarios where eccentricity effects could become observable and would give 
access to much interesting physics [37-42]. Third, eccentric inspirals in the extreme-mass- 
ratio regime will be key sources for a future mHz-band detector in space [43-47]. 

A gauge-invariant quantity for eccentric orbits, suitable for synergistic studies, was intro¬ 
duced by Barack and Sago in Ref. [24] (henceforth BS2011). This quantity is a straightfor¬ 
ward generalization of Detweiler’s redshift, obtained by averaging the time component of the 
particle’s four-velocity with respect to proper time over one epicyclic period of the motion. 
In other words, it is the ratio between the period measured in the coordinate time of a static 
observer at inhnity and the proper-time period. This “averaged redshift,” denoted here (U), 
is dehned with the dissipative piece of the GSF ignored. The functional relationship between 
(U) and the two invariant frequencies that characterize the motion is a gauge-invariant di¬ 
agnostic of the conservative eccentric-orbit dynamics. BS2011 calculated (U) (numerically) 
through linear order in the mass ratio for a sample of strong-held orbits, but they stopped 
short of attempting a calculation in a weaker-held regime where a meaningful comparison 
with PN results might be possible. The method of BS2011, which is based on a time-domain 
numerical integration of the relevant held equations, was best suited for tackling strong-held 
orbits, and its performance deteriorated fast with increasing orbital radius because of the 
longer evolution time required. 

Here we extend the range of BS2011’s calculation into the weaker-held regime, derive a 
3PN-accurate formula for (U), valid for any mass ratio, and compare between the numerical 
GSF results and the analytical PN prediction in the small mass-ratio limit. This is the hrst 
such comparison for noncircular orbits. It shows a good agreement for large and medium 
separations, and allows us to assess the performance of the PN expansion all the way down 
to the innermost stable orbit. Moreover, we are also able, through hts to the numerical GSF 
data, to extract some information about the 4PN approximation. 

Our numerical GSF calculation improves on that of BS2011 in both accuracy and weak- 
held reach. This improvement is achieved in two ways. First, our computation is based on 
the frequency-domain approach of Akcay et al. [48], in which the held equations are reduced 
to ordinary diherential equations. This ohers signihcant computational saving, particularly 
at lower eccentricities (e < 0.4). Second, we have found a way to signihcantly simplify the 
expression given in BS2011 for (U) as a function of the two orbital frequencies. The new 
form requires a simpler type of numerical input, which can be obtained at greater accuracy. 

This paper is organized as follows. In Sec. H we review relevant results for bound motion 
in Schwarzschild spacetime and for the redshift as dehned for circular orbits. We then extend 
the dehnition to eccentric orbits and obtain a simple expression for the generalized redshift 
(U) in terms of calculable perturbative quantities. Section HI discusses the numerics and 
sources of error, and displays a sample of numerical results for {U). In Sec. IV we perform a 
detailed derivation of the PN expression for (U) through 3PN order. Our calculations rely 
crucially on the known 3PN near-zone metric and the 3PN quasi-Keplerian representation 
of the motion. The numerical GSF and analytical PN results are compared in Sec. V. In 
Appendix A we establish the equivalence between our simplihed formulation of (U) and that 
of BS2011. Appendix B derives some useful PN formulas valid in the test-mass limit. 
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Table I summarizes some of our notation, for easy reference. In the GSF context, we 
denote the mass of the background Schwarzschild geometry by m 2 and the mass of the 
orbiting particle by mi, with the assumption that q = mi/m 2 ^ 1. In the PN context, the 
two particles of masses mi and m 2 have an arbitrary mass ratio q. We will set G = c = 1, 
except in Sec. IV where we keep these constants explicit in PN expressions. We use a metric 
signature (-, +, +, +). 


mi 

particle’s mass 

m 2 

black hole’s mass 

m = mi + m 2 

total mass 

q = mi/m 2 

mass ratio 

u = mim2/m? 

symmetric mass ratio 

A = {m 2 — mi)/m 

reduced mass difference 

nr 

radial (epicyclic) frequency 


average azimuthal frequency 


TABLE I: Important symbols. 


II. GENERALIZED REDSHIFT: FORMULATION IN SELF-FORCE APPROACH 


A. Bound geodesic orbits in Schwarzschild spacetime 


We hrst review relevant results for bound geodesic motion in Schwarzschild spacetime. 
We consider a test particle of mass mi moving on a bound (timelike) geodesic orbit in 
the Schwarzschild spacetime of a black hole of mass m 2 . Using Schwarzschild coordinates 
{t,r,9,(()}, we label the position of the particle by ^"(ro) = (tp(ro), rp(ro), 6 'p(ro), 0p(ro)), 
with four-velocity Uq = dXp/dro, where Tq is a proper-time parameter along the geodesic, 
and the label ‘0’ indicates normalization with respect to the background (Schwarzschild) 
metric ( 7 °^, i.e., Without loss of generality, we conhne the motion to lie in 

the equatorial plane, i.e., 6p = 7i/2, such that Uq = 0. We parameterize the geodesics by 
the two constants of motion: the specihc energy S = —Uot and specihc angular momentum 
C = uo<t„ where uoa = 5 '°/ 3 “o- 

The geodesic equation of motion is given by = 0, where V 0/3 is the covariant 

derivative compatible with the background metric For the above setup, this gives 


dtp 

dro 

d0p 

dro 



8 


/(^p) ’ 

C 


^2-Ueff(rp;£2), 


( 2 . 1 a) 

( 2 . 1 b) 


( 2 . 1 c) 


where /(r) = 1 — 2 m 2 /r, and 14 fr(r;£^) = /(r) (1 + C ?is an effective potential for the 
radial motion. Bound (eccentric) geodesics exist for 2a/2/3 < 8 <1 and C > 2\/3m2. For 
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a given bound geodesic, the radial distance r’p(ro) is confined to a finite range 2 m 2 < r^in < 
^"p(t'o) ^ ■'"max < oo witli T max denoting periastron and apastron radii, respectively. 
These two turning-point radii can be mapped bijectively to {£, £}. Thus the pair {rmin, r ^ax } 
can also parameterize the family of bound geodesics. 

Another such pair is given by the dimensionless “semi-latus rectum” p and the “eccen¬ 
tricity” e, defined by 


p = 


2 r^ 


m2 (r^ax 


^min 



^max ^min 

e =- 

^max T ^min 


These relations can be inverted to yield the Keplerian-like formulas 


pm2 

^max Z 5 

1 — e 


pm2 
1 + e ‘ 


( 2 . 2 ) 


(2.3) 


We can further express the specific energy and angular momentum in terms of p and e by 
solving the equation = V^g(r;£^) at r = {r min , T ma x}. Using Eqs. (2.3), this yields 
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{p — 2 — 2e){p — 2 + 2e) 
p(p — 3 — e^) 



£ 


pm2 

^p-3-e^ 


(2.4) 


Following Darwin [49] , we parameterize the radial motion using the “relativistic anomaly” 
X via 


^^p(x) 


pm2 

1 + e cos X ’ 


(2.5) 


where x = 0 and x = ^ correspond to periastron and apastron passages, respectively. Using 
Eq. (2.1c) with Eq. (2.5), we obtain 


drn 


m 2 p 


,3/2 


p — 3 — 


dx {1 + ecosxf \ P-^-‘^ecosx' 
which, with the help of Eqs. (2.1a), (2.1b) and (2.4), also gives 
dfn 


( 2 . 6 ) 


m 2 P 


{p — 2 — 2e){p — 2 + 2e) 


dx (p - 2 - 2 ecosx)(l + ecosx)^ V p- 6 - 2 ecosx 

d 0 p _ I V 


dx Y p — 6 — 2e cos x 

The functions ro(x )5 fp(x) and 0p(x) are all monotonically increasing along the 
radial periods in coordinate and proper times are calculated, respectively, via 


. 2 ^ dt, 


drn 


Tro= I ^dx, Tro= I ^dx, 

/o dx Jo dx 


(2.7a) 

(2.7b) 
orbit. The 

(2.8) 


and the accumulated azimuthal angle between successive periastron passages is 


— 


'0 


dx 


dx = 4 


P 

p — 6 + 2e 


ellipK ( 


4e 


\p — 6 -I- 2e 


(2.9) 
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Here ellipK(fc) = (1 — k sin^ 0) is the complete elliptic integral of the first kind 

and subscripts ‘0’ serve to distinguish the geodesic values T^o, %o and $0 from their corre¬ 
sponding GSF-perturbed quantities to be introduced below. For any (p, e) we have $0 > 27r, 
hence the periastron advances. 

We can now define the radial and (average) azimuthal frequencies via 

= ( 2 . 10 ) 

The pair provides a gauge-invariant parametrization of eccentric orbits. It should 

be noted, however, that the mapping between (p, e) and (12^, is not bijective: there exist 
(inhnitely many) pairs of physically distinct geodesics of different {p, e} values but the same 
set of frequencies. This degeneracy, first noted in BS2011, was thoroughly studied in [50]. 
The phenomenon is a feature of orbits very close to the innermost stable orbit. Since in this 
work we focus on less bound orbits (for the purpose of comparison with PN theory), the 
phenomenon of isofrequency pairing will not be relevant to us. 

In the parameter space of eccentric geodesics, stable orbits are located in the region given 
by p > 6 -|- 2e. The curve p = 6 -|- 2e is called the separatrix. Along it both $0 and Tro 
diverge, but Qfp remains hnite. This gives rise to the so-called “zoom-whirl” behavior [51], 
where the orbiting particle zooms in from far away, whirls around the black hole many times, 
thus accumulating a large azimuthal phase, then zooms back out. In the limit p —)■ 6 -|- 2e, 
the particle sits exactly at the peak of the effective potential and whirls inhnitely on an 
unstable circular geodesic. 


B. The redshift invariant for circular orbits 

Now let the particle’s mass mi be hnite but small, i.e.. 


q = mi/m 2 < 1, (2-11) 

and consider the ehect of self-interaction on the motion through 0{q). Within the context 
of linear perturbation theory, Detweiler and Whiting [52] showed that such a particle follows 
a geodesic motion in a certain smooth, ehective, locally-dehned spacetime with metric 

9al3 = 9% + ^a/3 • (2.12) 

Here, is a certain smooth piece of the physical (retarded) metric perturbation produced 
by the particle. The physical perturbation itself is a solution of the linearized Einstein 
equation, sourced by the particle’s energy-momentum, with suitable “retarded” boundary 
conditions. How may be computed in practice, on a Schwarzcshild background, is 
discussed, for example, in Ref. [23]. 

Within linear perturbation theory, may be split into a dissipative piece and a conser¬ 
vative (time-symmetric) piece, and the ehects of the two pieces may be considered separately. 
The conservative part of the perturbation is dehned as = 2 where 

and '' is a smooth perturbation constructed just like but starting with 

the particle’s “advanced” metric perturbation. Replacing —)■ in the effective 

metric (2.12) amounts to “turning off” the dissipation. The resulting equations of motion 
capture only conservative aspects of the dynamics. 
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In Ref. [15] Detweiler considered a particle in circular geodesic motion in the “conserva¬ 
tive” effective spacetime 

9T = 9% + hlp\ (2.13) 

In the absence of dissipation the orbit remains circular, and the spacetime possesses a helical 
Killing vector held, which, on the orbit, is proportional to the 4-velocity = da;“/dr. We 

introduce here r as a proper-time parameter along the geodesic in the effective metric 
with u°‘ normalized with respect to that metric, i.e., = —1. Thanks to the helical 

symmetry, all components of the particle’s 4-velocity are invariant under gauge transfor¬ 
mations that respect the helical symmetry [21]. Detweiler proposed to use the functional 
relationship between m* and D = as a gauge-invariant benchmark for the conservative 

self-force effect beyond the geodesic approximation. The frequency D is the circular-orbit 
reduction of the frequency dehned earlier for eccentric orbits. The quantity (or rather, 
its inverse) may be assigned a heuristic meaning of “redshift” (as measured in the smooth 
metric by a static asymptotic observer located along the helical symmetry axis), but 

it should be remembered that the true redshift, as measured in the physical metric of the 
particle, is, of course, divergent. 

Detweiler obtained [15] 

u\n) = m(,(D) + gMj,f(D), (2.14) 

where Mq = [l — 3(m2D)^'^^] is the geodesic limit, and 

(2.15) 

is the 0{q) correction arising from self-interaction. Note that the correction is dehned for 
a hxed value of D at the background, which ensures its gauge invariance. In Ref. [15] and sub¬ 
sequent work [16, 17] (see also [18-20]), Detweiler and collaborators calculated numerically 
the post-geodesic correction Mggf(D), and showed that it is consistent with corresponding PN 
expressions in an overlapping domain of validity. 

Detweiler’s numerical results were derived using the Regge-Wheeler gauge. An indepen¬ 
dent calculation using a direct numerical integration of the Lorenz-gauge form of the pertur¬ 
bation equations later recovered the same invariant relation M*gj-(D) [21]. This comparison 
highlighted a subtlety in the notion of invariance as applied to Mgg£(D): the gauge transfor¬ 
mation between the Lorenz-gauge metric perturbation and the Regge-Wheeler one does not 
leave Ugg£(D) invariant, due to a certain minor gauge irregularity of the Lorenz-gauge metric 
(that was hrst identihed in Ref. [53] and further discussed in [21]). Specihcally, the physi¬ 
cal metric perturbation does not vanish at inhnity when expressed in the Lorenz gauge; see 
Eq. (2.23) below. While the perturbation remains helically symmetric, the transformation to 
an “asymptotically flat” gauge like Regge-Wheeler’s (or the harmonic gauge of PN theory), 
in which Eq. (2.15) applies, has a generator that itself does not have a helical symmetry. As 
a result, the transformation introduces a correction to Mgg£(D). Denoting by the Lorenz- 
gauge metric perturbation, one hnds [21] 

• (2-16) 

The parameter a is extracted from the Lorenz-gauge perturbation as prescribed in Eq. (2.23) 
below; for a circular orbit it reads a = q{m 2 ^y^^UQ. One must be mindful, when working 
in the Lorenz gauge (as we do here), to take proper account of this gauge irregularity. We 
shall return to this point in more detail when discussing eccentric orbits. 


7 


C. The redshift invariant generalized to eccentric orbits 


Now consider an eccentric orbit snbject to the conservative effect of the GSF. In absence 
of dissipation, the orbit remains bonnd and has a constant radial period and a constant 
accumnlated azimnthal phase <I> per radial period. Hence it possesses a well dehned pair of 
freqnencies defined via Eq. (2.10) with the snbscripts ‘0’ dropped. The fnnctional 

relation between these invariant freqnencies and any gange-dependent set of parameters can 
be written as the snm of a “geodesic” term and a GSF correction; snch relations were derived 
in explicit form in BS2011 bnt will not be needed here. 

The GSF-pertnrbed orbit is a geodesic in the effective metric with 

tangent four-velocity m" normalized in It is easily checked that is no longer gauge- 

invariant in a pointwise sense when the orbit is noncircular. Instead, BS2011 suggested to 
consider the orbital average 

((7> ^ ^ Ft* dr = F , (2.17) 

'r Jo ' F 

where %. is the radial period measured in proper time r. BS2011 argued that (t/) is invariant 
under gauge transformations that respect the periodicity of the orbit and are well behaved 
(in a certain sense) at inhnity. We may split {U) in the form 

(G)(G,) = (G)o(G,) + q (G)gsf(a), (2.18) 

where Gj = {U)q is the geodesic limit of {U) taken with fixed 12*, and q{U)gsi is 

the GSF correction, defined for fixed 12*. The functional relation {17)gsf(12j) is an invariant 
measure of the GSF effect on the eccentric orbit, and it is the quantity that we will use for 
our GSF-PN comparison in this paper. 

The geodesic limit of (U) is given by 

= (2.19) 

IrO 

where the periods T^o and Tro may be calculated via (2.8) given the parameters p, e of the 
geodesic orbit. BS2011 describes a practical method for (numerically) inverting the relations 
12j(p, e) in order to obtain p(12j) and e(12j). This method may be used in conjunction with 
Eqs. (2.8) and (2.19) in order to compute {U)^ for given frequencies 12*. 

Our goal now is to express (17)gsf(12*) explicitly in terms of calculable perturbative quanti¬ 
ties (the metric perturbation and/or the GSF). Since hxing 12* Exes T*., the only contribution 
to (17)gsf (12*) comes from the 0{q) difference %—%o- From the normalizations g^/^UQU^ = —1 
and (p°^ -I- = — 1 one obtains 

^ = l + + (2-20) 

where terms of 0{q^) and higher are omitted. Since the contraction automatically picks 
out the conservative piece of the label ‘cons’ becomes redundant and we have dropped 
it. Neglecting subleading terms in the mass ratio g, we now obtain 

% - r.„ = (^1 - dr = ), 
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( 2 . 21 ) 


where the average is taken with respect to r (or tq) over a radial period. The 0{q) pertur¬ 
bation of {U) = T^-ITr at fixed hlj therefore reads 

9 {% - W = \(U),{ht)- (2.22) 

I'rOj ^ 

This would be our dual result for (17)gsf if were to be calculated in a suitable “asymp¬ 
totically flat” gauge. Our calculation, however, will be performed in the Lorenz gauge, which 
suffers from the aforementioned irregularity at inhnity. Let us now describe this irregularity 
more specihcally. For either circular or noncircular orbits, the Lorenz-gauge metric compo¬ 
nent htt tends to a hnite nonzero value at r —)• cx) (other components are regular). This 
behavior is due entirely to the static piece of the mass monopole perturbation, and there¬ 
fore the asymptotic value of hu does not depend on the angular direction even for eccentric 
orbits; it depends only on the orbital parameters. To remove this gauge artifact, following 
BS2011 we introduce the normalized time coordinate t = (1 -f a)t, where t denotes the 
original Lorenz-gauge time coordinate, and a = a{Qi) is given by 

a =-^htt{r ^ oo). (2.23) 

This normalization, which amounts to an 0{q) gauge transformation away from the Lorenz 
gauge, corrects the asymptotic behavior. Under t t we have, at leading order, 

hi ^ hi + 2ag°,{U)l = - 2a£(U)„. (2.24) 

Thus, to re-express {U)gsf in Eq. (2.22) in terms of the Lorenz-gauge perturbation, we need 
simply replace —)■ + 2aS{U)Q. We hnally get 

9 (U)^, = + a£{U)l . (2.25) 

Equation (2.25) is one of our main results, giving (U)gsf in terms of quantities directly cal¬ 
culable using existing GSF codes: the i?-field in the Lorenz gauge, and the corresponding 
asymptotic parameter a. It is clear that Eq. (2.25) reduces to (2.15) in the circular-orbit 
limit. As in the circular case, the expression for (U)gsf involves only the i?-£eld along the 
orbit (and the parameter a), and not the GSF itself. Our result (2.25) is much simpler than 
the one derived in BS2011 using a different procedure. In that work, certain simplihcations 
that reduce the expression for (f/)gsf to the form (2.25) have been overlooked. In Appendix 
A we establish the equivalence between the two results. 

III. NUMERICAL CALCULATION OF THE GENERALIZED REDSHIFT 

We have used the frequency-domain computational framework of Ref. [48] in order to 
compute (U)ggf for a large sample of orbits, focusing primarily on obtaining weak-held data 
for PN comparisons. Our calculation is based on Eq. (2.25), which takes as input the regular¬ 
ized Lorenz-gauge metric perturbation evaluated along the orbit (as well as the asymptotic 
value a, also to be read off the Lorenz-gauge perturbation). Since the GSF correction q{U)gsi 
(defined at hxed frequencies Gj) is of 0{q), it is sufficient to use as input the metric per¬ 
turbation calculated along geodesic orbits. For convenience we shall use p, e (as defined in 
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Sec. II A), rather than fij, to parameterize these geodesics, and will thus express our results 
in the form (t/)gsf = {U)gsi{p,e). It is important to emphasize that our results refer to the 
GSF correction to the functional relation (t/)gsf(f2j) defined for fixed invariant frequencies 
Qi, even though we use the geodesic parameters p and e as independent variables. These 
two facts should not be confused. 


A. Details of numerics and sources of error 

We use the eccentric-orbit GSF code of Ref. [48] to obtain the metric perturbation 
along the geodesic orbit. This code employs a frequency-domain approach, coupled with the 
method of extended homogeneous solutions of Ref. [54] , to compute the regularized metric 
perturbation It then outputs cit 2400 evenly spaced points along the orbit, and 

interpolates the numerical data using Mathematica's Interpolation function. In its default 
setting. Interpolation hts cubic polynomials between successive data points. Since 
is very smooth this level of interpolation is sufficient for our purposes. We subsequently 
calculate the orbital average using NIntegrate with the appropriate numerical inte¬ 

gration options/controls offered by Mathematica. The coefficient a is extracted, using (2.23), 
from the static monopole piece of the metric perturbation, whose construction is prescribed 
in App. B of [48]. Since this piece is essentially known analytically (its computation involves 
the evaluation of a certain orbital integral, easily done with Mathematica at extremely high 
accuracy), numerical error in our calculation of {U)gsi comes entirely from the numerical 
evaluation of h^„(x). Reference [48] contains a thorough analysis of error sources for h^„(x) 
and the GSF. Here, we briefly review two dominant sources. 

Each Fourier mode of our computation has associated with it a frequency, oj = 
where n and fh are integer harmonic numbers. The dominating source of numerical error 
depends on the value of u. For modes of sufficiently large frequency (m 2 a; > 10“^), the 
dominant error comes from the estimation of the contribution from the tail of uncomputed 
multipoles of large I values. Typically, we compute the contributions to all the /-modes 
up to and including / = 15, and estimate the remaining contribution to the mode sum by 
htting numerical data to suitable power-law models of the large-/ behavior [23, 48, 55]. This 
is a relatively well-modelled and well-controlled source of error, and it can be reduced in a 
straightforward manner using additional computational resources. 

For modes with small frequencies, m 2 a; < 10“^, a second source of numerical error takes 
over. This comes from rounding errors introduced when inverting the matrix of amplitude 
coefficients as part of the procedure for computing inhomogeneous solutions to the Lore n z- 
gauge held equations [48]. When u is very small, the matrix becomes nearly singular, and 
its inversion using machine-precision arithmetic introduces large errors. The problematic 
“nearly-static” modes occur generically in our calculation, since, given any orbital param¬ 
eters, there will exist values of m and n in the Fourier sum for which m 2 00 is very small. 
In practice, the sum over n and fh is truncated once our results reach a desired accuracy. 
Gonsequently, the problem is less severe for low-eccentricity orbits, where the effective fre¬ 
quency band is narrow, and more severe at high eccentricity, where the broad frequency 
band implies a higher chance of encountering nearly static modes. Ultimately, this restricts 
our calculation to orbits with eccentricities of e < 0.4. Low-ca modes are encountered also 
when the fundamental frequencies themselves are small, as with weak-held orbits—the main 
focus of the present work. Our code incorporates several methods for mitigating this small- 
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frequency problem (see Ref. [48] for details), but even with these techniques employed, our 
current calculation appears limited to orbits with p < 130; at larger p we observe a rapid 
reduction in accuracy. 

The issue of nearly-static modes has been addressed in a very recent paper by Osburn et 
al. [56], who proposed additional mitigation methods. These may be used to improve the 
performance of weak-held calculations in future work. 


B. Numerical results 

Table II displays a sample of our numerical results for {U)gsi- Parenthetical hgures indi¬ 
cate estimated error bars on the last displayed decimals; for instance, —0.0556761(1) stands 
for —0.0556761 ± 1 x 10“^. Additional data may be made available to interested readers 
upon request from the authors. Some of the data shown in the table are plotted in Fig. 1 of 
Sec. V, where we discuss the comparison with PN results. 

As a check of our frequency-domain computation, we compare our results for (f/)gsf with 
those obtained in BS2011 using a time-domain method. BS2011 provided a small sample of 
numerical results in the range p < 20 and e < 0.5. The comparison is shown in Table III. 
There is evidently a good agreement between the two sets of numerical results, although 
in some of the entries the values appear not fully consistent given the stated error bars 
(in all these cases the BS2011 values are smaller than ours). We have strong evidence to 
suggest that the source of disagreement is a slight underestimation of the magnitude of 
systematic error in the time-domain analysis of BS2011: We have tested the output of our 
frequency-domain code against accurate GSF data published in Ref. [56], and against yet 
unpublished redshift data calculated by van de Meent [57] (using a very different frequency- 
domain method based on a semi-analytical treatment of Teukolsky’s equation [58]). These 
comparisons strongly favor the frequency-domain data in the table. 

Also evident from the table is the fact that our code’s accuracy starts to degrade for 
e = 0.4; however, it still matches BS2011’s results to hve or six signihcant digits. No 
published numerical data exist to allow comparison beyond p = 20. (Reference [56] gives 
results for e < 0.7 and p < 90, but these are for the GSF components, not for (t/)gsf.) 


IV. GENERALIZED REDSHIFT: POST-NEWTONIAN CALCULATION 


We shall now derive the invariant relation within the PN approximation. 

Our calculations will be similar in spirit to those performed by Arun et al. [59, 60], except 
that we will consider the orbital average of a quantity that is related to the orbital dynamics 
of a binary of nonspinning compact objects, modelled as point particles, while Refs. [59, 60] 
calculated the orbital-averaged fluxes of energy and angular momentum radiated at inhnity. 
Furthermore, while these fluxes are invariant under the exchange 1 -f-)- 2 of the bodies’ labels, 
and requires knowledge of the gravitational held in the wave-zone, the generalized redshift 
(U) is a property of one particle, whose evaluation involves the near-zone metric. 
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e 

p = 10 

p = 15 

p = 20 

p = 25 

0.05 

-0.12878023(4) 

-0.07751154(5) 

-0.0556761252(8) 

-0.0434829334(1) 

0.10 

-0.1277540(3) 

-0.0768706(2) 

-0.05522166(7) 

-0.043132423(1) 

0.15 

-0.1260434(2) 

-0.07580395(6) 

-0.05446527(1) 

-0.042548963(6) 

0.20 

-0.123648(3) 

-0.07431376(9) 

-0.05340854(9) 

-0.041733648(3) 

0.25 

-0.120567(2) 

-0.07240329(5) 

-0.0520537(3) 

-0.04068802(6) 

0.30 

-0.1168020(6) 

-0.0700768(5) 

-0.05040377(4) 

-0.0394141(2) 

0.35 

-0.112352(2) 

-0.0673398(5) 

-0.0484623(4) 

-0.0379143(1) 

0.40 

-0.107221(2) 

-0.064199(1) 

-0.0462337(9) 

-0.0361916(2) 

e 

p = 30 

II 

CO 

II 

O 

o 

II 

0.05 

-0.0356833158(1) 

-0.0302606957(1) 

-0.0262706836(5) 

-0.0207905297(5) 

0.10 

-0.035398479(5) 

-0.0300209627(8) 

-0.0260637929(8) 

-0.0206282073(6) 

0.15 

-0.03492427(6) 

-0.029621799(3) 

-0.025719277(2) 

-0.0203578673(9) 

0.20 

-0.034261480(8) 

-0.029063791(5) 

-0.025237591(2) 

-0.019979802(2) 

0.25 

-0.03341121(1) 

-0.028347763(7) 

-0.024619373(3) 

-0.01949443(6) 

0.30 

-0.0323749(2) 

-0.0274748(4) 

-0.02386545(8) 

-0.01890227(5) 

0.35 

-0.0311543(2) 

-0.0264462(2) 

-0.0229768(4) 

-0.0182040(1) 

0.40 

-0.0297515(6) 

-0.0252634(2) 

-0.0219547(5) 

-0.0174004(4) 


p = 60 


p = 70 


p = 80 


p = 90 


0.05 -0.0172030750(2) 

0.10 -0.0170695612(3) 

0.15 -0.016847175(1) 

0.20 -0.016536122(4) 

0.25 -0.01613669(1) 

0.30 -0.01564925(2) 

0.35 -0.01507427(8) 

0.40 -0.0144123(3) 


-0.0146718447(3) 

-0.0145584700(2) 

-0.0143696130(3) 

-0.0141054255(4) 

-0.0137661204(6) 

-0.013351971(9) 

-0.0128633(1) 

-0.0123002(5) 


-0.0127901170(3) 

-0.0126916092(2) 

-0.0125275072(7) 

-0.0122979274(7) 

-0.012003033(2) 

-0.011643033(7) 

-0.01121819(3) 

-0.0107288(2) 


-0.0113362814(8) 

-0.0112491974(8) 

-0.0111041186(3) 

-0.0109011371(5) 

-0.010640382(2) 

-0.010322020(3) 

-0.00994625(3) 

-0.0095133(2) 


p = 100 


p = 110 


p = 120 


p = 130 


0.05 

0.10 

0.15 

0.20 

0.25 

0.30 

0.35 

0.40 


-0.0101792669(2) 

-0.0101012344(10) 

-0.0099712296(7) 

-0.0097893274(4) 

-0.0095556341(8) 

-0.009270281(4) 

-0.00893344(2) 

-0.0085453(2) 


-0.0092365822(2) 

-0.0091658975(1) 

-0.0090481316(3) 

-0.0088833462(1) 

-0.0086716263(10) 

-0.008413085(2) 

-0.00810786(1) 

-0.0077561(2) 


-0.0084537150(1) 

-0.0083891149(2) 

-0.0082814820(1) 

-0.0081308692(1) 

-0.0079373488(7) 

-0.007701017(5) 

-0.00742198(3) 

-0.0071004(2) 


-0.0077931956(9) 

-0.0077337155(1) 

-0.0076346120(1) 

-0.0074959280(1) 

-0.0073177296(5) 

-0.007100089(2) 

-0.00684311(3) 

-0.0065469(2) 


TABLE II: Numerical data for the GSF contribution (f7)gsf to the generalized redshift (defined 
with fixed invariant frequencies Oj), for various eccentric geodesic orbits in a Schwarzschild back¬ 
ground. The orbital parameters e (eccentricity) and p (semi-latus rectum) are defined in Sec. II A. 
Parenthetical figures indicate estimated error bars on the last displayed decimals. 
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e 

p = 10 

p = 15 

p = 20 

0.1 Here 

BS2011 

-0.1277540(3) 

-0.1277554(7) 

-0.0768706(2) 

-0.0768709(1) 

-0.05522166(7) 

-0.05522177(4) 

0.2 

-0.123648(3) 

-0.1236493(7) 

-0.07431376(9) 

-0.0743140(1) 

-0.05340854(9) 

-0.05340866(4) 

0.3 

-0.1168020(6) 

-0.1168034(6) 

-0.0700768(5) 

-0.0700771(1) 

-0.05040377(4) 

-0.05040388(4) 

0.4 

-0.107221(2) 

-0.1072221(5) 

-0.064199(1) 

-0.0641991(1) 

-0.0462337(9) 

-0.04623383(4) 


TABLE III: Our frequency-domain numerical results for (?7)gsf and the corresponding time-domain 
values from BS2011. Each cell shows our result (top) in comparison to BS2011’s (bottom). The 
relative disagreement between the two data sets is ~ 10“®, roughly consistent with the magnitude 
of error bars. As discussed in the text, evidence suggests that our results are accurate to within the 
error bars given, whereas the magnitude of error in the time-domain data is slightly underestimated 
in some cases. No time-domain data exist for p > 20 to allow comparison in the weak-field domain. 


A. Redshift variable in standard harmonic coordinates 

1. The regularized 3PN metric 


Throughout Sec. IV we assume that mi < m 2 , and treat mi as the “particle” orbiting the 
“black hole” of mass m 2 . The redshift of the particle can be computed from the knowledge 
of the regularized PN metric gagivi) = Qagii) Yi), generated by the two bodies and evaluated 
at the coordinate location yi(t) of the particle, as [16] 


U = u\ 



v?v 


1 



(4.1) 


where vf = (c, vi), with vi = dyi/dt the coordinate velocity of the particle. The generalized 
redshift will be given by the proper-time average of Eq. (4.1) over one radial period. 

The regularized PN metric gagivi) was itself computed up to 2.5PN order, in harmonic 
coordinates, in Ref. [61]. This calculation was then extended to 3PN order in Ref. [16], partly 
based on existing computations of the 3PN equations of motion using Hadamard regular¬ 
ization [62] and dimensional regularization [63]. Reference [16] performed two calculations 
of the 3PN regularized metric, using both Hadamard and dimensional regularizations, ob¬ 
taining the same metric but expressed in two different harmonic coordinate systems. The 
two metrics were found to differ by an inhnitesimal 3PN coordinate transformation in the 
“bulk,” i.e., outside the particle’s worldlines, and also by an intrinsic shift of these worldlines. 
Combining Eqs. (4.2) and (A15) of Ref. [16], the 3PN-accurate expression of the regularized 
metric reads, in the standard harmonic coordinates corresponding to the use of Hadamard 
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^ As usual we denote by ri 2 = |yi — y 2 | the coordinate separation, by ni 2 = (yi — y 2 )/ri 2 the unit direction 
from particle 2 to particle 1, and by V 12 = Vi — V 2 the relative velocity, where = dyo/dt is the 3-velocity 
of particle a. The Euclidean scalar product between two 3-vectors A and B is denoted {AB). Parentheses 
around indices are used to indicate symmetrization, i.e., A^'^bA = + A^B'‘). 
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(4.2c) 


Since we are considering only the conservative part of the binary dynamics, we did not in- 
clnde in (4.2) the dissipative 2.5PN radiation-reaction terms; these can be fonnd in Eqs. (7.6) 
of Ref. [61]. Notice the occnrence at 3PN order of some logarithmic terms, containing two 
constants r[ and (one for each body) that have the dimension of a length. These nltraviolet 
(UV) regnlarization parameters come from regularizing the self-field of point particles using 
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the Hadamard regularization of Ref. [62], The constants r'^ and r '2 are gauge-dependent, as 
they can be arbitrarily changed by a coordinate transformation of the bulk metric [16, 62] 
or by some shifts of the worldlines of the particles [63]. The metric coefficient 5 'oo( 2 /i) also 
involves a constant vq that originates from the infrared (IR) regularization of the metric at 
spatial infinity, as discussed in Ref. [16]. This arbitrary IR scale should also disappear from 
the final gauge-invariant results. 

We introduce the expression (4.2) of the regularized 3PN metric into the definition (4.1) 
of the redshift, and expand in powers of 1/c, keeping all terms up to 0(c~^). This gives an 
expression for U as a function of the two masses mi and m 2 , the coordinate separation ri 2 , 
and the scalar products (ui 2 Ui), (^ 12 ^ 2 ), (uiUi), (^ 1 ^ 2 ) and {V 2 V 2 ), as well as the regularization 
constants vq, r[ and r^, in an arbitrary reference frame. The resulting expression is too 
lengthy to be displayed here. 


2. Reduction to the center-of-mass frame 


We wish to specialize the previous expression to the center-of-mass (CM) frame, which 
is consistently defined at 3PN order by the vanishing of the center-of-mass integral deduced 
from the 3PN binary equations of motion [64]. This condition yields expressions for the 
individual positions and velocities relatively to the CM in terms of the relative position 
y = yi — y 2 and relative velocity v = vi — V 2 [65]. Since these results play an important role 
in our algebraic manipulations, we recall here the expressions for the functional relationships 
ya[y, v] in the harmonic gauge that was used to derive the regularized metric (4.2). Thus, 

yi = [X 2 + iy {Xi -X 2 )V]y + iy (W - X 2 ) Q v + 0 ( 0 -"*), (4.3a) 

y 2 = [-Xi + u{Xi- X 2 ) P] y + z/ (Xi - X 2 ) Q V + 0 ( 0 "®), (4.3b) 


where u = mim 2 /m^ is the symmetric mass ratio and Xa = ma/m, with m = mi -|- m 2 the 
total mass of the binary. The coefficients V and Q depend on the parameters m and u, the 
separation r = |y|, the relative velocity squared = (uu), and the radial velocity r = {nv). 
They explicitly read [65] 
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Again, we did not inclnde the radiation-reaction contribntions at 2.5PN order. A logarithmic 
term contribntes at 3PN order in (4.4a); it involves a particnlar combination Tq of the gange 
constants r[ and that is dehned by 

(Ai - X 2 ) In r' = In r[ - X^ In . (4.5) 


By compnting the time derivatives of Eqs. (4.3)-(4.4) and by applying, where necessary, 
an iterative order-reduction of all accelerations by means of the CM equations of motion 
given in Eqs. (3.9)-(3.10) of Ref. [65], we obtain expressions analogous to (4.3)-(4.4) for the 
particle’s individual velocities as functions of the relative variables y and v. Replacing the 
positions and velocities by their CM expressions ya[y,v] and Va[y,v] yields 3PN-accurate 
expressions for the scalar products (ni 2 ni), (^ 12 ^ 2 ), (uiUi), {V 1 V 2 ), (V 2 V 2 ) as functions of 
r, r and Finally, the CM expression for the redshift f/[r, r,n^] in standard harmonic 
coordinates takes the form 


1 + ^Cn+^ f^lPN + ^ f^2PN + ^ f^3PN + o(c ®) , 


where the various PN contributions read 


(4.6) 
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the test-mass limit of particle 1 corresponds to z/ —?• 0. Since the redshift (4.1) is a property 
of particle 1, the expressions (4.7) are not symmetric by exchange 1 -f-)- 2 of the bodies’ labels. 
The redshift of particle 2 is simply obtained by setting A —)■ —A in Eqs. (4.7). As expected, 
the regularization constants ro, r[ and that enter the expression (4.2) of the regularized 
3PN metric appear in the CM expression (4.6)-(4.7) for the redshift. In Sec. IV D we will 
check that the orbital averaging cancels out the dependance on these arbitrary length scales. 


B. Redshift variable in alternative coordinates 

In the previous section we obtained an expression for the redshift variable in the standard 
harmonic (SH) coordinate system, namely the coordinate system in which the 3PN equations 
of motion were originally derived [62, 65]. These coordinates are such that the equations of 
motion involve some gauge-dependent logarithmic terms at 3PN order. Importantly, these 
logarithms prevent the use of the 3PN quasi-Keplerian representation of the binary motion 
(reviewed in Sec. IV C below), thus impeding the averaging of the redshift over an orbit. 
Therefore, it is useful to have the expression for the redshift in a modihed harmonic (MH) 
coordinate system, without logarithmic terms in the equations of motion, such as the one 
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used in Refs. [59, 60]. Alternatively, we shall use ADM-type coordinates, which are also free 
of such logarithms at 3PN order in the equations of motion. Both the MH coordinates and 
the ADM coordinates are suitable for a 3PN quasi-Keplerian parametrization of the motion 
[66]. This will require us to re-express the redshift in terms of the variables r, r and in 
these alternative coordinate systems. 


1. Modified harmonic coordinates 


The trajectories y'ait) of the particles in MH coordinates are related to their counterparts 
jait) in SH coordinates by some 3PN shifts ^a{t) of the worldlines induced by a coordinate 
transformation in the “bulk,” namely y'^ = ya+^a [62]. Therefore, in the CM frame, the MH 
coordinate separation y' is related to the SH coordinate separation y through y' = y + ^, 
where the relative shift ^ — ^2 is given by [59] 


^(SH MH) 


22 , 
y c6r2 ^ 



n + o(c '^), 


(4.8) 


with n = y/r the unit direction pointing from particle 2 to particle 1. Following [59], we 
introduced the “logarithmic barycenter” Tq of the constants r[ and rg, (not to be confused 
with the IR constant ro): 

In Tq = Xi In -|-X 2 In r 2 . (4.9) 

The expression U'[r, r, n^] for the redshift in MH coordinates can then be deduced from the 
formula for U[r, r, n^] in SH coordinates by means of the functional equality U' = U + 6^U, 
where 

, (4.10) 

with 


= «) + (P(e"), 


= {ni) + 




-«) + ci(e^), 


= 2(y) + (P(e^). 


(4.11a) 

(4.11b) 

(4.11c) 


Since the relative shift (4.8) comes at 3PN order, the nonlinear terms in Eqs. (4.10) and 
(4.11) contribute at leading 6PN order, and can thus be neglected. Plugging the expression 
(4.8) into Eqs. (4.11), we hnd the explicit expressions 


ikkm 
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MH) 


22 ( r 

In — 




SH -S. MH) 


22 G^m^v 
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44 G^mfiv 
3 


2 r In ( — 


- 3r2) In 


(4.12a) 

(4.12b) 

(4.12c) 


In order to compute the change 5^U in the redshift induced by the relative shift (4.8), we 
only require the Newtonian expression for f/[r, r, n^], which is given by Eq. (4.7a). Combined 
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with (4.10) and (4.12), this gives 
11 


(Wk 




(1 + A — 2z/) + ^(1 + A) (v^ — 3r^) — 

+ ^(1 + A — 3z/ — A n) — 3r^) — (1 + A — 2z/)^ In | + o(c“®). (4.13) 


Adding the above shift to the formnla (4.6)-(4.7) for the redshift in SH coordinates yields the 
expression for the redshift in MH coordinates. Since 17' = t/ + is a fnnctional eqnality, 
the resnlting MH redshift is expressed as a fnnction of the “dnmmy” variables r, r and . 

Adding together Eqs. (4.7d) and (4.13), we hnd that the UV regularization constant 
disappears from the expression for 17'[r, r,n^] in MH coordinates. However the UV and IR 
constants r( and tq remain and enter the result through the logarithmic contributions 



11 G^rrfiv / 2 


37^ 


y){»-—-'"(a 



(4.14) 


For circular orbits, such that r = 0 and = Gm/r + 0{c~‘^), these logarithmic contributions 
cancel out. We will see that the factor (n^ —3r^ —Gm/r)/r^ vanishes when averaged over one 
radial period, such that the constants ro and r( will cancel out from the hnal, gauge-invariant 
result for the orbital-averaged redshift, as expected. 


2. ADM-type coordinates 


Similarly, the individual trajectories y(j(t) of the particles in ADM coordinates are related 
to the trajectories ya{t) in SH coordinates by some shifts ^a{t) of the worldlines: y'^ = ya+^a 
[64, 67]. In the CM frame, the ADM coordinate separation y' is related to the SH coordinate 
separation y through y' = y + ^, where the relative shift ^ — ^2 reads [59, 65] 
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from which the authors of Ref. [59] deduced, using Eqs. (4.11), the transformation of variables 
that we need to compute the redshift in ADM coordinates:^ 
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^ The remainder in Eqs. (4.10) and (4.11) is of order 4PN, which is still negligible in the transformation 
to ADM coordinates. 
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The expression G'[r, r,n^] for the redshift in ADM coordinates can then be deduced from 
the result (4.6)-(4.7) in SH coordinates via the functional equality [/' = [/ + 6^U. Using 
the expressions (4.10) and (4.16), the SH redshift is found to be modihed by 2PN and 3PN 
corrections that read 
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(4.17) 


Although the additional contribution (4.17) in ADM coordinates is more involved than its 
counterpart (4.13) in MH coordinates, they share the same logarithmic terms. Thus, adding 
Eqs. (4.7d) and (4.17) we find that the UV regularization constant disappears from the 
expression for f/'[r, r,v^] in ADM coordinates, while the constants rg and r( remain and 
enter the final result through the logarithmic terms (4.14). 


C. The generalized quasi-Keplerian representation 

Before we discuss the orbital averaging of the redshift in Sec. IV D, we must summarize 
the 3PN generalized quasi-Keplerian (QK) representation of the motion of Memmesheimer 
et al. [66]. Indeed, since averaging over one radial period is most conveniently performed 
using an explicit solution of the equations of motion, the generalized QK representation is 
an essential input for our 3PN calculation. The QK representation was originally introduced 
by Damour and Deruelle [68] to account for the leading-order IPN general relativistic effects 
in the timing formula of the Hulse-Taylor binary pulsar. It was later extended at 2PN order 
in Refs. [69-71], in ADM coordinates, and more recently at 3PN order [66] in both ADM 
and harmonic coordinates. 

We first introduce the mean anomaly 


i = Qr(t- fper) , (4-18) 

where fper is the coordinate time at a periastron passage and Qr = 2ti/T r is the radial fre¬ 
quency (also known as the mean motion n), i.e., the frequency associated with the periodicity 
Tr of the radial motion. The mean anomaly simply maps one radial period t G [fper, ^per + Tr) 
to the trigonometric interval £ G [0,27r). We then adopt a parametric description of the bi¬ 
nary’s motion in polar coordinates, in the CM frame, in terms of the eccentric anomaly 
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u G [0,27r). At 3PN order, this parametrization reads 

r(M) = (1 — cosm) , (4.19a) 

£{u) = u — CiSinti + ft sin id + gt{V — u) + itSin2V + htsinSid, (4.19b) 

0('u) = 0per + K {y + f^sm 2V + (70 sin ?>V + sin 4 V" + sin bV) , (4.19c) 

where 0per is the value of the orbital phase when t = tper, at a periastron passage, K = 1 + k 
is the fractional angle of advance of the periastron per orbital revolution, such that the 
angle of return to periastron is given by 4) = 27tK (equivalent to A<l> = 27rk), and the true 
anomaly V is dehned by 


V{u) = u + 2 arctan ( — ^ ^ , (4.20) 

\ 1 — P 0 cos u J 

with y = [l — (1 — /e 0 . Equations (4.18)-(4.20) provide a 3PN-accurate general¬ 

ization of the usual Keplerian representation of the Newtonian motion.^ 

The previous generalized QK representation is complete only once the orbital elements 
fir, K, ttr, Ct, Cr, 60 , ft, Qt, it, ht, f^, Q^, and ^0 have been related to the hrst integrals of 
the motion, namely the binding energy E and the orbital angular momentum J, both per 
reduced mass jjt = mim 2 /m. Following Ref. [59], we shall instead make use of the convenient, 
dimensionless, coordinate-invariant quantities'^ 


2E . _ 2EJ^ 

^ {Gmy 


(4.21) 


such that 6 > 0 and j > 0 for a generic bound eccentric orbit (since E < 0 for such orbits). 
Notice the PN scalings 6 = (P(c“^) and j = (P(c°). Therefore, we shall consider expansions in 
powers of the PN parameter e, with coefficients depending on j and u. In ADM coordinates, 
the 3PN-accurate expressions for the orbital elements read [59, 66 ] 
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(4.22b) 


^ In the Newtonian limit, is the semi-major axis, the three eccentricities coincide (e* = Cr = = e), an 

eccentric orbit does not precess [K = 1), and ft = gt = H = ht = f 4 > = 94 , = 14 , = h^j, = 0. 

^ For circular orbits, we have the well-known Newtonian limits e ~ jt? ^ Gm/(r(?) and j ^ 1. 
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The eccentricities and are all such that j = 1 — + 0{c ^) at Newtonian 

order; they start differing from each other at leading IPN order. 

The expressions (4.22) are specihc to the ADM coordinates. Before we give the corre¬ 
sponding expressions in MH coordinates, let us recall an important point related to the use 
of gauge-invariant variables. As shown in Ref. [69], the functional forms of D,. = 2ti/T r and 
K = ^/ (27r) as functions of gauge-invariant variables like £ and j are identical in different 
coordinate systems. In particular we have the exact same relations in MH coordinates as in 
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ADM coordinates: 


^MH _ ^ADM ^ ^ ^4_23a) 

j^MK ^ ^ADM ^ ^ _ (4.23b) 


We may therefore use any combination of and K instead of the constants of the motion 
£ and j to parameterize in a physically meaningful way a given eccentric orbit (assuming a 
one-to-one relation). Following Ref. [59], we introduce the frequency which is 

a natural generalization of the circular-orbit frequency D,® and we dehne the dimensionless 
coordinate-invariant parameters (remember that k = K — 1) 

V c3 ) ’ T 



The PN parameter x is 0{c~‘^), while i is merely Newtonian at leading order (the relativistic 
periastron advance hrst appears at IPN order). The choice of variables (4.24) is the obvious 
generalization of the gauge-invariant variable x that is commonly used for circular orbits. It 
will thus facilitate checking the circular-orbit limit. In Sec. IV E, we shall express our hnal 
results in terms of either of the two sets of gauge-invariant parameters {e,j) or {x, i). 

To compute the invariant relationships {U){e,j) and {U){x, l) from the expressions (4.6)- 
(4.7) and (4.10) for the redshift in MH coordinates, we shall also need expressions for the 
orbital elements ar, Ct, e^, ft, gt, it, ht, g^, i,/, and in these coordinates. They are 
given by Eqs. (4.22c)-(4.221), to which we must add the differences [66] 
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® Note that coincides with the average angular frequency of the motion: = (f) = f- f(t) dt. 
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Notice, in agreement with the comment made earlier in Sec. IV B 2, that the MH coordinates 
differ from the ADM coordinates at leading 2PN order. 
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D. Orbital average of the redshift 


We are finally in a position to compnte the generalized redshift 


1 


cTr 


{U) = Y U(T)dT, 

>r Jo 


(4.26) 


which coincides with the ratio T^jTr of the coordinate time period and the proper time 
period %. of the radial motion.® The averaged redshift (4.26) can be written in the convenient 
alternative forms^ 

1 If 1 r2TT ^ p2-K pi / \ 

‘ ^ (4.27) 


= 

u can be c 
averaging in MH coordinates. 


Tr Jo U{t) 271 Jo U{i) 277 Jo U(u) 

where f = d^/dn can be computed from Eqs. (4.19b) and (4.20). We first perform the orbit 


1. Orbital average in MH coordinates 

Using the generalized QK representation (4.18)-(4.20), (4.22)-(4.23) and (4.25), the vari¬ 
ables r, r and = r^-|-r^0^ that enter the expression (4.6)-(4.7) and (4.13) for the redshift in 
MH coordinates can be expressed as functions of the binding energy e, the time eccentricity 
Ct = and the eccentric anomaly u. The integrand in Eq. (4.27) then reads 
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E tX 

(1^ 

vr=-l ^ 


aNieus) , o / Mn(l-etcosn) 

-^ + V ^Niet, e) - -— . 

et cos uy^ (1 — e* cosn)^'^ 


(4.28) 


The computation of the coefficients avr and is straightforward, but the resulting expres¬ 
sions are too cumbersome to be reported here. The integral in (4.27) is readily performed 
thanks to the following formulas, which are valid for all integers N ^ 1 \ 
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y + ^/y"^ 


(4.29a) 

(4.29b) 

(4.30) 


We note that the logarithmic contributions in (4.28) arise at 3PN order from those terms 
proportional to In (r/rg) and In {r/r'J) in Eq. (4.14). Indeed, combining Eqs. (4.19a), (4.22d), 
(4.22e), (4.25b) and (4.25c), one finds 


In ( - 1 = In 

To 


— ) In (1 — Cf cosu) -|- 0(c ^), 
ro. 


(4.31a) 


® Beware that, although we are using the same symbol to denote the generalized redshift in Eqs. (2.17) and 
(4.26), the former definition is restricted to linear order in the mass ratio, while the latter holds for any q. 
^ Notice the simple relation {U)^{l/U)^ = 1, where (•),. (resp. (•)() denotes an averaging over one radial 
period with respect to the proper time r (resp. the coordinate time t). 
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(4.31b) 



Hence some coefficients ajv in (4.28) depend on the regularization constants tq and r[ through 
In [ttr/roj and In . However, when averaged over one radial period, these terms cancel 

out from the hnal expression, because they appear only through the vanishing combination 


2/2(e) - 5 / 3 ( 0 ) + 3(1 - e^)h{e) = 0 . (4.32) 

The hnal expression for {U){e, e^^) is thus free of the regularization constants tq and r(. 

Implementing all the above integrations, the expression (4.6)-(4.7) and (4.13) for the 
redshift in MH coordinates can be averaged over an orbit. Up to 3PN order, the generalized 
redshift (4.26) then takes the form 

(U) = 1 + e + , (4.33) 


where the PN coefficients depend on the symmetric mass ratio z/, the reduced mass difference 
A = \/l — 4z/, and the time eccentricity et in MH coordinates (hence et = e^^). They read 
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(4.34d) 


For notational simplicity we did not add a label on et to indicate that it is the time eccen¬ 
tricity in MH coordinates. (No such label is required over e, which is gauge invariant.) This 
point should be remembered when comparing expressions derived in different gauges, as we 
shall do next. 
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2. Orbital average in ADM eoordinates 


We shall now perform an independent calculation in ADM coordinates. We start from the 
expression for the redshift in ADM coordinates, as given by (4.6)-(4.7) and (4.17), employ 
the appropriate QK parametrization and perform the orbital averaging as outlined above. 
We hnd that the form (4.28) is obtained also in the ADM case, with the same coefficients 
Pn but different coefficients in general. The result for the generalized redshift in ADM 
coordinates is of the form 


(t;) = 1 + + M2™ + M3™" £" + o(£‘), 


(4.35) 


where the various coefficients depend on u, A, and the time eccentricity in ADM coordinates 
(hence et = ef-^^), and read 
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(4.36d) 


Although the coefficients (4.34) and (4.36) coincide through 2PN order, the 3PN coefficients 
and are different. A useful internal check of the PN calculations of the generalized 

redshift in MH and ADM coordinates is the verihcation that the equality of Eqs. (4.33)- 
(4.34) and (4.35)-(4.36) holds if and only if the time eccentricities and are related 
by 

1 + 17.^ d + 0(£=') 
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(4.37) 
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This relation is in perfect agreement with what is predicted from using different QK repre¬ 
sentations of the motion, namely Eq. (4.22d) together with (4.25b). 


E. Gauge-invariant formulations 


To compare the analytical PN predictions with the numerical results of the GSF calcula¬ 
tion (Sec. Ill), it is best to use a coordinate-invariant relationship. We shall thus replace the 
coordinate-dependant time eccentricity Ct in favor of the coordinate-invariant angular mo¬ 
mentum variable j. Substituting the PN expansion (4.22d) into Eq. (4.36), or alternatively 
Eqs. (4.22d) and (4.25b) into (4.34), we get 

{U) = 1 -|- Wn £ + WiPN + l^2PN + WsPN + o{£^) , (4.38) 


where 
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(4.39a) 

(4.39b) 
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Since the relationship {U){e,j) is coordinate-invariant, it is physically meaningful. How¬ 
ever, the binding energy E and angular momentum J are not easily accessible to perturbative 
GSF calculations, so a direct comparison is not obvious. Thanksfully, Eq. (4.38) can also 
be expressed using the invariant parameters (4.24) defined with respect to the fundamental 
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frequencies and Indeed, inverting the PN expansions (4.22a) and (4.22b) yields 
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We thus have the leading-order relationships a; = £ + (P(c ^) and t = j-|-C>(c ^). Introducing 
the expansions (4.40) into Eq. (4.38)-(4.39), our main PN result reads 
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where the various PN coefficients, which depend on the variable l as well as on the particle’s 
masses, read up to 3PN order 
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The noncircular nature of the motion only explicitly enters the result at leading IPN order 
via the invariant parameter l. Since we have the qualitative behavior l ~ 1 —e^, this suggests 
that the effect of the eccentricity on {U) will be moderate (at least in the weak-held regime). 


1. Circular-orbit limit 

Another key check of the results (4.39) and (4.42) is provided by the circular-orbit limit. 
For such orbits, the two constants of the motion are no longer independent variables. Indeed, 
the angular momentum variable, say j©, is related to the energy e by the 3PN gauge-invariant 
expansion [72] 
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(4.43) 


It can be checked that the eccentricities e*, e^, all vanish when j is replaced by (4.43) in 
Eqs. (4.22d)-(4.22f) and (4.25b)-(4.25d). The invariant result (4.38)-(4.39) then reduces to 
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Setting Ct —)■ 0 in Eq. (4.34) or (4.36) yields the same expression. 

We then replace the constant of the motion £ in favor of the frequency-related parameter 
X [recall Eq. (4.24)], using the well-known 3PN-accurate expression for the binding energy 
as a function of the circular-orbit frequency, namely [see, e.g., Eq. (232) of Ref. [14]] 
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Finally, replacing e in (4.44) using (4.45), we recover the known 3PN result for the circular- 
orbit redshift (see Eq. (4.10) of Ref. [16]): 
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(4.46) 


Interestingly, at Newtonian order, the averaged redshift (JJ) along an eccentric orbit has the 
same functional form as Uq in the case of a circular orbit. This shows that the effect of the 
eccentricity cancels out at Newtonian order, because of the orbital averaging. 

Alternatively, we can also combine Eqs. (4.40b), (4.43), (4.45) to obtain the PN expansion 
of the invariant relation iQ{x) in the circular-orbit limit, namely 
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(4.47) 


and introduce this expression into (4.41)-(4.42) to recover (4.46). 

Our third and last check of the correctness of the formula (4.42) will be to recover the 
known result in the test-particle limit. 


2. Extreme mass-ratio limit 


The 3PN result (4.41)-(4.42) is valid for any mass ratio q = mi/m 2 . To extract from this 
result the contribution due to the conservative piece of the GSF, we introduce an alternative 
set of dimensionless coordinate-invariant parameters, better suited than (x, l) to the extreme 
mass-ratio limit g 1: 
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We substitute the relations x = g)^/^ and z, = A (1 -|- g)^/^ in (4.41)-(4.42), and expand 

in powers of the mass ratio g, neglecting terms of 0{q^) or higher. The 3PN result for the 
sum of the test mass, GSF and post-GSF contributions reads 
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In the test-particle limit g = 0, we recover the 3PN expansion of the fully relativistic result 
(2.19) for a geodesic orbit, as derived in App. B. The 3PN prediction (4.50c) could be 
compared with future calculations of the second-order GSF [73-78]. 

Finally, we may express the result (4.50b) for the 3PN expansion of the GSF contribution 
to the generalized redshift by means of the usual parametrization of bound timelike geodesic 
orbits in Schwarzschild in terms of the semi-latus rectum p and eccentricity e (see Sec. II A). 
Substituting for y and A from Eqs. (Bl) into (4.50b), we find 
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where je = 1 — e‘^. For small eccentricities, we may write 


(G)gsf = a-|-&e^-|-ce^-|-(ie®-|- 0{e ^), 


(4.52) 


where the weak-field expansions of the coefficients a{p), b{p), c{p) and d{p) read 
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and higher-order terms in the eccentricity all contribute at leading 2PN order. 


(4.53a) 

(4.53b) 

(4.53c) 

(4.53d) 


V. COMPARISON OF POST-NEWTONIAN AND SELF-FORCE RESULTS 

In Fig. 1 we plot our data for (U)gsf as a function of p for a sample e = {0.1, 0.2, 0.3, 0.4} 
of eccentricities. We show, superposed, the corresponding IPN, 2PN and 3PN predictions 
from Eq. (4.51). The insets display the relative differences between the GSF data and the 
successive PN approximations. We make the following observations: 

(i) There is an excellent agreement between the numerical GSF results and the analytical 
PN prediction at “large” p, in what should be considered a very strong test of both 
calculations. This is a first demonstration of such an agreement for noncircular orbits. 

(ii) The PN series appears to converge uniformly to the GSF result at any p for any fixed 
e in our survey, at least through 3PN order. 
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(iii) The 3PN formula reproduces the GSF results extremely well even in what might be 
considered a “strong-held” regime: at p = 10 it does so to within ~ 1% for e = 0.1 
and to within a few percent for e = 0.4; at p = 20 the agreement is already at the 
level of one part in a thousand. 




(a) e = 0.1 


(b) e = 0.2 



FIG. 1: Numerical GSF output for (17)gsf (black data points) versus analytical PN approximations 
(solid curves). Each panel shows (P)gsf as a function of semi-latus rectum p for a fixed eccentricity 
e. Insets display, on a log-log scale, the relative differences A”™ = |1 —t/nPN/(P)gsf|) where C/nPN is 
the PN approximation through nPN order. In both the main plots and the insets, the three curves 
correspond, top to bottom, to the IPN, 2PN and 3PN approximations. Solid curves in the insets 
are the analytical PN residues 1 — Pipn/Pspn (upper curve) and 1 — t/ 2 PN/P 3 PN (middle curve); 
for the lower curve we have fitted the simple model 1 — P 3 PN/(P)gsf = P”"^ («! -b 0:2 Inp + as/p). 


We can make the comparison more quantitative by attempting to extract the large-p 
behavior of the numerically computed function (17)gsf(p, e). Our strategy will be to fit the 
numerical data against the PN model (4.51), leaving the numerical coefficients as unknown 
htting parameters, later to be compared with the analytically known values. Given the 
relative sparseness of data available, we shall not attempt a simultaneous £t over p and e, 
but rather fit over each of the two dimensions separately, as described below. We will follow 
a “marginalization” procedure, whereby each of the PN orders is htted for in turn, assuming 
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the analytic values of all terms at lower PN order. Since the circular limit of (t/)gsf has been 
computed previously at great accuracy [17, 18, 22, 34], we are able to accurately “remove” 
the circular (e-independent) part of (t/)gsf from the data, htting only for the e-dependent 
residue. This should allow to £t the eccentricity-related terms of interest here with greater 
accuracy. 

Let us now describe this procedure in more detail. We assume the e-expanded form (4.52) 
of the full PN expression (4.51). The term a{p) is the circular-orbit limit of (f/)gsf, which 
has been computed to at least ten signihcant hgures in Refs. [17, 18, 22, 34]. By subtracting 
off these numerical data from ours, we construct a new data set for the difference 

= {U)gsi - a{p) = b{p) + c{p) + • ■ ■ . (5.1) 

We assume that the functions b{p),c{p),... admit expansions in p~^ as in Eqs. (4.53), but 
pretend that the PN coefficients are unknown: 

b = p~^ + bi p~‘^ + &2 P~^ H-, 

c = cip~'^ + C 2 P~^-\ -, (5.2) 

where subscripts are mnemonics for the PN order at which coefficients occur, and we have 
hxed the “Newtonian,” 1/p term of b{p) at its known value of unity. Our goal is to determine 
the coefficient bn, Cn, ... from the numerical data for To this end, we hrst prepare 

subsets of data where in each subset p is hxed and e varies. We £t each subset with respect 
to e using the model (5.1), including terms through This yields three one-dimensional 

data sets, representing b{p), c{p) and d{p). 

Focusing hrst on the data set for b{p), we ht it against the PN model 

N 

b{p) = p~^ + {bn + bn^lnp) , (5.3) 

n=l 

in which = 62 °® = 63 °® = 0 , since logarithmic terms are known not to occur before the 
4PN order [79, 80].® The truncation order N is left as a control parameter; by varying it we 
obtain a rough estimate of the numerical uncertainty in the htted values of the parameters. 
We apply a marginalization procedure, whereby to determine we set all bn'<n at their 
known analytic values. We use this procedure to estimate the values of 61 , 62 and 63 , and 
we later similarly determine Ci. Our results are shown in Table IV, alongside the known 
analytic values for these parameters. We see a good agreement through 3PN order in the 
0{e^) term, and at IPN order in the 0{e^) term. 

Unfortunately, the accuracy of our current code (and its limited utility at e > 0.4) does 
not seem to allow us an accurate extraction of &n> 4 , Cn> 2 , or any of the The reason for 

this can be appreciated from Fig. 2, where we compare the amplitudes of the 3PN and 4PN 
terms with the amplitude of numerical noise in our (U)gg| data. Note that, while the “signal” 
from the 63 term lies well above the noise, the C 2 signal is buried deep inside it. Since our 
data is limited to relatively small eccentricities, it is clear why we have less “handle” on the 
Cn [0{e^)] terms than on the bn [0{e^)\ terms. 

® The form of the circular-orbit limit, in which the PN expansion of (C/)gsf is known analytically up to a 
very high order [17-20, 81], suggests that the function b{p) could also involve powers of Inp. However, 
those would contribute at even higher orders than we consider here, so we do not include them in the PN 
model (5.3). 
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Coefficient 

Estimate 

Exact result 

61 

-^4.0002(8) 

-f4 

b2 

-h7.02(3) 

+7 

bs 

-14.5(4) 

-14.312... 

Cl 

- 2 . 00 ( 1 ) 

-2 


TABLE IV: Best-fit values for the PN coefficients and Cn [Eqs. (5.1) and (5.2)] as extracted from 
the numerical data, compared to their known exact values. Parenthetical figures are estimated 
fitting uncertainties in the last displayed decimals, obtained by varying the value of the truncation 
index N in the htting model [e.g., Eq. (5.3) for b{p)]. The exact value of 63 is — (5/3 + 4l7r^/32). 



EIG. 2: Absolute magnitude of various PN terms (“signal”) compared to the magnitude of numer¬ 
ical error in the data (“noise”), shown as a function of p for e = 0.1 (left panel) and e = 0.3 (right 
panel). The red (upper solid) curve shows the 3PN term of the O(e^) piece of {U)gsi, and the 
blue (lower solid) curve shows the 2PN term of the O(e^) piece. Comparison with the magnitude 
of numerical noise (black dots) suggests that 63 should be easily discernible while C2 might not. 
This is confirmed by attempting to fit the data against PN models, as detailed in the text. The 
dashed curve estimates the amplitude of the 4PN term of the O(e^) piece of (V)gsf, which is not 
known analytically. We used here the values 64 = —1500 and 6^® = 250 chosen from the middle 
of the estimated range shown in Eqs. (5.4). This 4PN signal appears to lie just over the noise and 
is detectable. However, as it can be seen from the near overlap of the densely dashed (green) and 
sparsely dashed (brown) curves, the I 64 I and 64 °® terms become almost equal in magnitude as p 
increases, hence making it very difficult to extract the individual values of 64 and 64 °®. 


Figure 2 also suggests that we might have just enough “signal” coming from the 0{e‘^) 
terms at 4PN to allow a rough estimation of the coefficients 64 and 64 °^, which are not known 
analytically. We have experimented htting to a large number of models of the form (5.3), 
where all the analytically known coefficients are pre-specihed, and varying both the cutoff N 
and the number of nonzero logarithmic terms. We hnd that htting uncertainties are almost 
as large as the htted values themselves. However, we are able to conhdently constrain the 
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values of 64 and 64 ® to lie within the ranges 


-2000 <64 <-1000, 

+ 150<6^<+350. 


(5.4a) 

(5.4b) 


Future analytic calculations of the 4PN terms may be checked against these predictions. 

Our current code does not allow the determination of unknown PN coefficients related 
to eccentricity with any greater accuracy. To improve on our predictions would require (i) 
to push the reach of the computation to higher eccentricities and larger p, and at the same 
time (ii) to reduce the numerical error in the calculation of (tl)gsf- Some improvement may 
be achieved using the method of Ref. [56] , which is a slightly more advanced variant of our 
method. More significant improvements may have to await the development of eccentric- 
orbit GSF codes based on the Teukolsky equation [58, 82]. We expect such codes to start 
delivering accurate numerical results in the very near future. 
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Appendix A: Equivalence between our expression for (U)gsf and BS2011’s 

BS2011 give their formula for (U)gsf in their Eq. (84). Adjusting notation and rearranging 
the terms, their expression reads 



where 



(A 2 ) 


and AX denotes the GSF correction to a quantity X, holding hxed p and e, rather than the 
invariant frequencies Dj. BS2011 give explicit expressions for AT^, A$ and A% in terms of 
GSF quantities, but these will not be needed here. We observe that the expression (Al) is 
much more complicated than our result, Eq. (2.25). Our goal here is to show that the two 
expressions are, in fact, identical. 
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(A3a) 

(A3b) 


To this end, we note the two key relations 

C* = £ {U)l , 

a = ^ [(s - c%) {u}„ - 1], 

which shall be derived below. Snbstituting these Cr and Cip into Eq. (Al) and using Tro = 
TroKU)^, we obtain 

5W„(/d) + o£{(/>^ + -^(£Ar,-£A4.-Ar.). (A4) 

For this to be identical to Eq. (2.25), the sum of three term in brackets on the right-hand side 
should vanish. Indeed, writing (dr/dy)^ = — 5 fQ,/ 3 (dx"/dx)(dx^/dx) and perturbing linearly 
with A, holding p, e and Vp (or, equivalently, p, e and y) hxed, we hnd 

A(dr/dx) = £ A{dt/dx) - C A{d(j)/dx) , (A5) 

which, upon integrating over a radial period, gives 


ATr = £ATr-CAd>. (A6) 

Hence Eq. (A4) reduces to our Eq. (2.25) for (t/)gsf- 

It remains to establish the relations (A3a) and (A3b). This can be achieved by manipu¬ 
lating the explicit elliptic-integral representations of fir and {U)p, given in BS2011, but 
this approach involves much ungainly algebra and will not be presented here. A much neater 
derivation uses general results derived from the Hamiltonian formulation of geodesic motion 
in Kerr spacetime [83]. Start by averaging UQUoa = — 1 with respect to t over a radial period 
of the geodesic orbit, to obtain 

{U)~^ =£- - fir Jr , (A7) 

where Jr = (27r)~^ § uor dr = (27rTo)~^ ('^^o)^ is invariant action variable (per mass 

mi) associated with the radial motion [84]. This relation is the Schwarzschild reduction of 
Eq. (3.4) of Ref. [83]. In addition, we require a relation between the partial derivatives of £, 
C and Jr with respect to flj. The necessary relation follows most directly from the general 
variational formula (“hrst law”) 


5£ = VLs5C + VLr6Jr 


(AS) 


established in [83] [this form is the reduction of Eq. (3.5) therein to Schwarzschild spacetime, 
with a hxed black-hole mass m 2 , and with suitable notational adjustments]. Here 5£, 6C 
and 6 Jr correspond to an arbitrary variation of a geodesic with frequencies flj onto a nearby 
geodesic. If we regard £, C and Jr as functions of flj, we obtain 


d£ Q o 

dVLi ^ dfli '' dfli 


0 . 


(A9) 


Taking the partial derivative of (A7) with respect to fl^ and using (A9) immediately leads 
to (A3a). Equation (A3b), in turn, is obtained by taking the derivative of (A7) with respect 
to fir, then using Eq. (A9), and hnally substituting for Jr from (A7). 

The above establishes the equivalence of our simple expression (2.25) and the BS2011 
result (Al). The simplihcation obtained here owes itself primarily to the two key relations 
(A3a) and (A3b), which have unfortunately gone unnoticed (by two of us) in BS2011. 
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Appendix B: Post-Newtonian expansion of {U)q 


Here we consider a test mass on a bonnd geodesic orbit aronnd a nonspinning black hole 
of mass m 2 and obtain the PN expansion of the relationship This calcnlation 

provides a powerfnl check of onr PN result (4.41)-(4.42), because it is based on a different 
formalism, it makes use of an alternative parametrization of the motion, and it is performed 
using a different coordinate system. 

Since the relationships (2.7)-(2.10) cannot be inverted analytically to yield the expressions 
for the parameters p and e as functions of the frequencies 12 ^ = 27r/Tro and = $o/Tro, 
we shall work perturbatively, expanding all quantities in powers of the small parameter 1/p. 
From Tro and $0 we define the invariant parameters y= and X= 3y (4>o/27r — 1)”^ 

[recall Eq. (4.48)]. Expanding the formulas (2.7)~(2.10) up to 3PN order, we obtain 


y = -\i + 

p [ 

133 

IT ~ 


2(1-Je) 


P 


+ 




' je] + 5je 


VTe 


pz 


48je - 35jf' + 27j/ + 25jf/' - y + o{p-^) | , 


(Bla) 


A = je 


1 - 
1849 


11 7 , \ 1 

T + yd p+ 


75 23 , ^ , 0/2 73 .2 


pz 


849 , 


45 , 


17 , 


_i3/2 _ —I'S/a 

Je ' A Je 


95 


16 


2341 

192 


Je] ^ + o{p n 


po 


(Bib) 


where we introduced the notation if, = l — e^. In the limit of vanishing eccentricity, e —?• 0, 
we have the simple relation y = p~^ + o{p~‘^). Actually, we know that for circular orbits 
the relation y = 1/p holds exactly, such that in Schwarzschild coordinates the semi-major 
axis coincides with an invariant measure of the orbital radius. [This, however, is no longer 
true at 0{q) in the GSF approximation.] For circular orbits, the 3PN-accurate relationship 
between the invariants y and A then reads 


^ 1 9 9 2 27 3 , 3 , 

= 1 - 2^ “ 4^ “ + "(y) ■ 


(B2) 


Inverting the relations (Bl) yields expressions for the semi-latus rectum p and eccentricity 
e (or equivalently je = 1 — e^) as functions of the invariant parameters y and A. Up to 3PN 
order, we find 


1 

p 


Je 




A 


65 5 

^^47x 


4A, 

25 1 

^ ^ 4A2 


16 16A 

2255 A 
) 


8A2 J 

+ o{y^) 


,'7 llA / 5 63 


2 

16A2 ) ^ 


5 , 145 , 221 ^ 95 

6 ^ ^ m ^ 8A3/2 


289 


+ 


263 




64A2 192A3 J 


y^ + o{y^) 


(B3a) 


(B3b) 


We now have all the pieces required to compute the relation {U)Q{y, A) up to the required 
PN order. The generalized redshift is defined as 

1 rTrO / 1 /•27r I- ■■ \ -1 


(^)o = r- 

' r 


rO Jo 


( 1 

Uq{to) dro = — 

V-f rO 


dto dy 


dy 4(x) 


(B4) 
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where Tro is the proper time period of the radial motion. From the expressions (2.4), (2.5), 
(2.7a) and (2.8), we find 


P 




21 


(^)o “ ^ + 77 1 O ) 7 + ( ~ 6j, 






P 






249 

— Vje - 24je 


105j, 




(B6) 


Finally, substituting for (p, e) in terms of (y. A) in Eq. (B5), using (B3), we obtain the 
3PN-accurate coordinate-invariant relation 


6 


{u)o = i + ^y+ o + ^-Th/+-TF + 


2^ Vs A 

605 ^ 117 ^ 411 

^ 64v^ ^ ^ 


41 57 


16 


1755 21 

+ TTX + 


4v^ 

1797 


3 


37 
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15 A 
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4A3/2 2AV 


y 


20 


32A3/2 4A2 




In the circular-orbit limit, we may introduce the PN expansion (B2) for X{y) in (B6), expand 
in powers of y up to the appropriate PN order, and recover the 3PN expansion of the fully 
relativistic result Uq = (1 — 3y)~^^‘^. Although the result (B6) can in principle be extended 
up to an arbitrarily high PN order, we only need here the 3PN approximation to the exact 
result. Comparing with the formula (4.50a) derived from our 3PN calculation valid for any 
mass ratio, we hud perfect agreement. 
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